Multiscale science is to study the phenomena which couple disparate spatial or temporal scales.^[1] In the realm of computational material science, the properties of material systems should be considered on a large scale comprising 10 to 12 orders of magnitude, ranging from 1 nm to a few hundred meters.^[2] For different length scales (electronic, atomistic, microscopic, mesoscopic, and macroscopic scales), a great variety of well-established theories (quantum mechanics, molecular mechanics, continuum mechanics, classic mechanics, etc) can be provided as invaluable tools to treat the relevant phenomena on each special hierarchy. Nonetheless, almost every approach used on the relative scales has its own unavoidable limitations. Quantum mechanical (QM) calculations are indispensable for treating chemical reactions, charge transfer, electron excitation, and magnetism in materials,^[3] which, theoretically, could describe the behaviors of the material system on all the scales. However, mathematical complexity of this method is so great that it is often hindered by its significant computational cost, with most current applications being limited to less than 1000 atoms, and the limitations on model size impede the representation of many interesting geometries and can cause inaccuracies in energy and force due to mistreatment of the effect of long-range interactions of either electrostatic or elastic fields.

A less computationally demanding approach than the QM calculations, for example, is to use the molecular mechanics, implemented by the molecular dynamics (MD) simulations based on the empirical interatomic potentials. They can treat systems with millions of atoms and have been very successful in realistic applications. However, in many cases, they are inadequate, because they neglect the effect of electrons and thus lack crucial information about how the microstructure influences the macroscale behavior of the system. More importantly, the construction of a high-quality potential used in it is obtained by fitting the experimental or computational bulk data such as equilibrium volume, bulk modulus, elastic constants, cohesive energy, etc., which is very time consuming and unable to capture all of the complexity that is found during electronic QM calculations.^[4]

The constraints and inadequacy of the traditional monoscale approaches are precisely the motivation for the establishing of multiscale modeling. Multiscale algorithms hybriding several computational models specialized on disparate scales can offer a promising solution to achieving small-length-accuracy while involving the environmental effects on a much larger scale. They can distribute the computational power where it is most needed, and thus can make a reasonable balance between systematic accuracy and computational efficiency. However, the reasons for coupling different length scales are not only out of considerations for the computational resources. From the viewpoint of physics, many problems are inherently multiscale, that is, the processes on different scales are ruled by different physical laws, and several scales interact strongly to produce the observed macroscale behaviors of the materials. In order to find out what happens exactly on each scale, we should use the right tool for the right part of the system.^[5] Thus, the essence of multiscale modeling is to develop better physics models. Multiscale approach offers a more unified view to modeling, and opens up unprecedented opportunities.

Based on the level of how different scales interact with each other, the multiscale approaches can be divided into two categories: sequential and concurrent methods.^[6] When the interaction is comparatively weak, sequential multiscale algorithm is used, which is the one that the consuming-computational-power on a finer-length scale is employed to obtain parameters for use in a more approximate or phenomenological methodology on a coarser-length scale.^[7–9] Sequential coupling has been constrained mainly to the circumstances when the needed information is just a small number of parameters or functions of very few variables. For this reason, sequential approach is often referred to as parameter passing. One landmark work is that in the 1980’s, Clementi et al.^[10] used sequential multiscale method to predict tidal circulation in Buzzard’s Bay. In the 1990’s, a sequential multiscale modeling strategy was used by Wang et al.^[7,11] to study the macroscopic properties such as relation between impurity and stacking-fault energy from the electronic structure of an impurity–stacking-fault complex by using the tight-binding Hamiltonian and recursion method. Another example is that Wang et al.^[12] employed the transition state theory combined with the QM approach to calculate the migration and vacancy formation energies of Ni and Al atoms in the pure and Re-doped Ni-based single crystal superalloys. The corresponding results were then transmitted to the kinetic Monte Carlo method to estimate the atomic diffusivity values of Ni and Al atoms under different temperature conditions.

In contrast to the sequential paradigm, the concurrent method is necessary for the system whose behavior on one scale depends strongly on what happens to the others. The concurrent approaches tend to couple algorithms of disparate scale hierarchies simultaneously with hybriding procedures performed in some overlapping domain in order to resolve the important multiscale physics. Typically, concurrent methods combine two distinct length scales together, or three at most.^[6]

During the decades, various concurrent multiscale methods have been developed. Generally they can be categorized into energy-based and force-based approaches.^[13] Energy-based coupling strategy is to develop a well-defined global energy functional from which the forces are calculated.^[9] For covalently bonded organic molecules, Karplus^[14] and others were awarded the Nobel prize in chemistry 2013 for the development of hybrid models that combine the advantages of classical and quantum methods to describe complex chemical systems. Energetic coupling term that models the interaction between the classical and the quantum systems was constructed. For covalently-bonded bulk materials, Broughton et al.^[5] proposed a multiscale algorithm coupling continuum (finite element [FE] method) to statistical (molecular dynamics [MD] method) to quantum mechanics (tight binding [TB] method) to examine crack propagation in silicon. A Hamiltonian was defined for the entire system as H_tot = H_FE + H_FE/MD + H_MD + H_MD/TB + H_TB. “MD/TB” and “FE/MD” indicate handshake regions between different methods. As for the MD/TB region, dangling bonds of atoms of the TB region were saturated with univalent atoms. For the FE/MD boundary, FE and MD approaches each contributed half the weight to the handshake Hamiltonian. Following a trajectory dictated by this Hamiltonian would result in a conserved total energy and numerical stability. For simple metallic materials, Choly et al.^[15] proposed a multiscale strategy coupling density functional theory-based quantum simulations to classical atomistic simulations with the expression of the energy of the whole system as E[I + II] = E₁[I] + E₂[II] + E_int[I, II], where E₁[I] was calculated with DFT and E₂[II] was estimated with a classical potential. The crux lies in the concurrent coupling E_int, which is treated quantum mechanically via the orbital-free density functional theory. In this case, the degree of freedom is the electron charge density in the DFT region, and the total energy functional is directly minimized with respect to the charge density.

Force-based coupling strategies use different energy functionals and forces for different domains.^[9] The total energy is not clearly defined in this method, but it is defined by starting from a definition of the coupling in terms of a physically-motivated prescription of the forces on unequilibrated atoms.^[13] Woodward et al.^[16–19] proposed the first-principles lattice Green function boundary condition method (FP-GFBC) coupling the strain field produced by a single dislocation to the long-range elastic field with using a flexible boundary condition. In this method, the incompatibility forces produced within the supercell containing the dislocation were adapted by relaxing all the atoms of the system according to the Greens function solution. Up to now, this method has mainly been used in the metallic system. Based on the FP-GFBC approach, Liu and Wang^[20] calculated the solute–dislocation interaction energy in and around the dislocation core and then predicted the strengthening of single crystal superalloys by adopting an analytic model for the strength, or stress to move a dislocation, owing to the random field of solutes.^[21] Nair et al.^[4] developed a concurrent multiscale method by coupling a region governed by quantum mechanics to a surrounding region governed by continuum mechanics. This method is based on the coupled atomistic discrete dislocation (CADD) framework by Shilkrot et al.^[22] with the atomistic region replaced by a quantum-mechanical atomistic region. So it is referred to as QM-CADD, which can be applied to non-metals and metals.

The disadvantage of the energy-based strategies is that it is difficult to eliminate the non-physical ghost forces^[23] at the boundary between different domains. On the other hand, the force-based coupling schemes can be slow to equilibrate, converge to non-equilibrium states, and become non-conservative. The review paper^[24] has pointed out that the present force-based approaches have superior accuracy over the energy-based approaches. However, it seems unlikely that the energy conservation could easily be restored into the force-based scheme. Therefore, there is a necessity of developing an energy-based scheme that has buffer regions and we need to find a proper method to partition the energy spatially in a physically meaningful way.

In this paper, we report an energy-based concurrent multiscale approach, the energy density method (EDM)^[25–27] and its recent development by our research group. We show that the EDM can provide a reasonable solution to the above problems. First, the EDM has employed a transition region to smoothly connect the separated QM and MD domains. Second, the EDM could partition the QM and MD energies spatially to each atomic site, and we refer to this quantity as the site energy. Third, we introduce the weighting parameters into the transition region, representing the contribution of the QM and MD site energy to each atomic site. Based on the above consideration, we can give a unified expression of the energy functional of the entire system that combines the QM and MD descriptions. Moreover, the EDM is used to study the multilayer relaxation of the Ni (001) surface structure and is validated against the periodic DFT calculations. The results show that the concurrent EDM could combine the accuracy of the DFT description with the low computational cost of the MD simulation and is suitable to the study of the local defects with the influence of the long-range environment.

The basic idea of the EDM approach is to spatially divide the whole system into three distinct regions denoted as a QM, a MD, and a transition (TR) region (see Fig. 1(a)), respectively. The QM region includes vacancies, impurities, dislocations, surfaces, or other defects that may cause substantial electronic redistribution, and we use a DFT method to treat this area. The MD domain represents the long-range electrostatic or elastic field, providing an outer environment for the QM region. The MD region is close to the perfect crystal lattice, and thus the accuracy of a MD method would be sufficient. Since we are most interested in the metallic systems, and the embedded-atom method (EAM)^[28] has proven to be a significant improvement in simplified total-energy calculations for metals, we choose the EAM-based MD method to deal with the MD region. A third region spatially located between the above two domains hybridizes the different physical laws of two different scales, and hence is called the transition region.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. (color online) (a) Side view and (b) top view of atomic structure of Ni (001) surface in a 10 × 10 × 30 supercell with 12000 atoms; (c) schematic diagram of the side view of Ni (001) surface. Layer number increases for deeper layers from surface to bulk. The Ni atoms are indicated by blue balls.

The energy of the total system is explicitly written as a sum,

In the following, we introduce the concept of the site energy into the DFT method and the MD method respectively, which could meet the above three requirements. In the DFT method, the total energy can be separated into band structure energy or one-electron energy E_band (sum of the values of one-electron eigenvalue ), and the so-called electrostatic term E_es,^[29] i.e.,

After the EAM was brought up, Daw^[32] derived a simplified expression for the cohesive energy of the metallic system of EAM from approximations to density functional theory. Thus the quantum mechanical background of the EAM and the close relationship between the two are proven. The expression of the cohesive energy of the DFT method

Based on the above discussion, equations (7) and (8) can be decomposed to each atom as follows:

We assume that the cohesive energy of the QM and TR region obtained by EDM should be equal to that calculated by using the DFT approach alone. We specify this proposition as the energy constraint condition

Next, ionic forces are obtained by differentiating the total energy. Giving these forces, we could perform the ionic relaxations and the system can be driven into a stable equilibrium structure by minimizing the ionic forces on all the ions in the system. We choose the powerful molecular dynamics leap-frog algorithm. The following algorithm is iterated:

.

Fcc Ni is the basic matrix of a high-temperature Ni-based superalloy.^[33] Here we use the EDM to study the surface structure of the fcc metal Ni. Since the EDM is suitable for the defect systems with long-range environment, a simple case is to perform the surface calculation. It is necessary to have a microscopic understanding of the surface relaxation among the basic questions of surface science.^[34] First-principles calculations of multilayer relaxations are still computationally expensive due to the large size of the unit cell necessary to simulate the surface structure. The potential used in the EAM is obtained by fitting the experimental data such as equilibrium volume, bulk modulus, elastic constants, cohesive energy, etc, and hence, the EAM studies are not so reliable as the first-principles calculations.^[35] The multiscale paradigm EDM could provide a useful algorithm in dealing with such problems. In the following, we illustrate how the EDM is performed in the multilayer relaxation of Ni (001) surface.

The EDM is performed in a 10 × 10 × 30 supercell along the x [100] × y [010] × z[001] direction. The periodic boundary conditions are implemented in the x and y direction as shown in Fig. 1. The atoms in each (002) plane having the same z coordinate are equivalent in the surface structure relaxation, and therefore share the same site energy and weighting parameter in the energy density method.

For the QM and TR region, the first principles calculations are performed within the generalized-gradient approximation of Perdew–Burke–Ernzerhoff^[36] for the exchange–correlation functional as implemented in the Vienna ab initio simulation package.^[37] The DVM-DAC software^[38] is utilized to calculate the electrostatic term of the quantum mechanical site energy . The cut-off energy of atomic wave function is set to be 350.5 eV for Ni. After considering the symmetry of the atoms on this perfect surface (without vacancy or other defects), the QM and TR regions are set to have 7 and 18 (002) layers with 14 and 36 atoms along the z direction, respectively, and a special k-point sampling method is utilized for the integration by setting 25 × 25 × 1 with Monkhorst–Pack scheme^[39] in the first irreducible Brillouin zone. Convergence criterion employed for the electronic self-consistency is set to be 10⁻⁵ eV. The ionic relaxation is stopped when the forces on all of the atoms are less than 0.02 eV/Å. The EAM potentials of Ni developed by Du et al.^[40] are adopted for the MD calculations.

The site energies are obtained during the EDM procedure, and the results reveal that could be seen as an energetic quantity representing the characteristic of the lattice structure. The weighting parameters ω_l for each (002) atomic layer in the TR region are also obtained. From boundary QM/TR to TR/MD, gradually decays from 1 to 0, while increases from 0 to 1. We perform the full relaxation of all atomic positions. The relaxation procedure should be achieved when the maximum force is less than 0.02 eV/Å. After 26 iterations steps, the force convergence criterion is reached as shown in Fig. 2.

Fig. 2.

	Figure Option ViewDownloadNew Window
	Fig. 2. (color online) Variations of maximum factor with the number of iteration step, showing convergence of maximum force on the atoms in the entire system in each iteration by using EDM and FLWMS to simulate the Ni (001) surface.

For comparison, we perform the ionic relaxation by using the force linear weighted mixing scheme (FLWMS),^[9] that is, the contribution of the forces calculated by the DFT method to the atomic force in the TR region is weighted according to a linear function that varies from 1 at the interface QM/TR to 0 at the interface TR/MD, while by MD, from 0 to 1. The FLWMS is a typical force-based concurrent multiscale method. From Fig. 2, we can see that the convergence process of EDM is more stable and gentle than that of FLWMS. For the FLWMS, it will take about 32 iteration steps before the maximum force is less than 0.02 eV/Å. The convergence speed of the EDM is almost 6 iteration steps (∼ 19%) faster than that of the FLWMS. Note that if nonlinear effects like complicated structural and chemical defects are taken into consideration, the disparity of the results between the two methods would be acute. Thus, our EDM has a more intrinsic nature to deal with the multiscale problems.

In order to assess the performance of the EDM calculation, we also relax the QM and TR region together by using the DFT method, and the entire system by using the MD method for comparison. The results are shown in Fig. 3.

Fig. 3.

Figure Option ViewDownloadNew Window

Fig. 3. (color online) Calculated multilayer relaxation Δdn,n + 1 = 100(dn,n + 1 − d0)/d0 of Ni (002) surface by using EDM. The relative surface atomic structure is given in Fig. 1. DFT calculational models for QM and TR regions are both surface slab consisting of 25 (002) atomic layers repeated periodically along the z direction, separated by 12 Å vacuum layers. Periodic boundary condition is implemented in x and y directions.

In Fig. 3, d_n,n+1 is the interlayer spacing between the atomic (002) layers n and n + 1 (n = 1, 2, 3,..., 49) along the z direction as illustrated in Fig. 1(c). So d₁₂ is the difference between the z coordinations of atoms in the first and second (002) layers. The result of the surface multilayer relaxation Δd_{n,n + 1} indicates the change in the interlayer spacing between layers n and n + 1 as compared with the bulk interlayer spacing given in percentage Δd_{n,n + 1} = 100(d_{n,n + 1} − d₀)/d₀, where d₀ = a₀/2 = 1.76 Å. Positive sign indicates the outward expansion while negative value corresponds to the inward contraction.

In the QM domain and at QM/TR interface (layers 1–7), because of the quantum effect of the surface structure, the DFT calculations could give a reliable estimation for the relaxation of this region. The most important concern in this paper is the quantum effect and relaxation results of the relatively narrow area of a few atomic layers near the surface. For complex surface defects, the DFT method could not give correct answers by using a small model comprising scant top most atomic layers near the surface. The empirical interatomic potential that has not specifically performed the surface properties fittings is not suitable to treating such problems. When other surface defects, namely, adatoms, vacancies, steps, etc. are concerned, it will be much more complex for the empirical interatomic potentials to handle. The results shown in Fig. 3 suggest that our EDM model could give the values for the surface area (layers 1–7) which are comparable to the DFT calculations and thus has the ability to maintain the DFT accuracy in the key surface region while considering the further bulk environment. Moreover, oscillatory multilayer relaxation is observed in both DFT and EDM results in Fig. 3. In the bulk Ni, the ion cores are screened by the conduction electrons around them. The creation of a surface changes the coordination number and electron density of atoms in the first layer. This makes the distribution of electrons turn into corrugation at the surface.^[41] The electrons will spread themselves out, leading to a smooth charge density at the surface. This generates the asymmetric screening of the ion cores in the first layer and a net electrostatic force which pushes them into the crystal, and thus reducing d₁₂. This phenomenon is called the Smoluchowski effect.^[42] Shortening d₁₂ stabilizes atoms in the first layer and destabilizes atoms in the second layer. The atoms in the second layer can be stabilized by moving the atoms in the third layer away, that is, by expanding d₂₃. The oscillatory multilayer relaxation results can be viewed as the simultaneous response following the alteration of d₁₂.^[43]

In the region around TR/MD interface (layers 18–25), it turns already into bulk Ni without the influence of the surface. So there should be little variation of Δd_{n,n + 1}. The DFT calculations applied to the QM and TR regions could not give the right result due to the factitious interception of the bulk Ni. While the efficient MD method applied to the entire system has no such questions. In Fig. 3, we can see that the EDM model could obtain a similar atomic structure to the MD result in this region, which means that the EDM could effectively reduce the boundary effects between adjacent domains and to some extent has achieved the seamless coupling.

For comparison, in Fig. 4, we also present the results of the surface multilayer simulated by FLWMS. Through the analysis of Fig. 4, we can see that the EDM and FLWMS calculations have a similar degree of accuracy, while the energy-based EDM simulation has no problems such as the non-conservative energy mentioned above.

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) Comparison between calculated multilayer relaxation Δdn,n + 1 values of the Ni (002) surface obtained by using EDM and FLWMS. The relative surface atomic structure is given in Fig. 1.

We have reported a new concurrent multiscale approach, the energy density method, which couples the quantum mechanical and the molecular dynamics simulations simultaneously. The coupling crossing different scales is achieved by introducing a transition region between the QM and MD domains. In order to construct the energy formalism of the entire system, the concept of site energy and weight parameters of disparate scales are discussed. Moreover, the EDM is applied to the study of the multilayer relaxation of the Ni (001) surface structure. Our results reveal that the concurrent EDM could successfully combine the accuracy of the DFT description with the low computational cost of the MD simulation to deal with surface problems with long-range environment. The EDM could also effectively reduce the boundary effects between adjacent domains, and to some extent has achieved the seamless coupling. The energy-based EDM has a similar degree of accuracy to FLWMS, while it does not have problems such as the non-conservative energy. It should be further employed to explore the complicated systems involving dislocations or cracks, in which there are couplings between the local defects and the long-range strain fields.